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SUMMARY 


Approximate nonlinear inviscid theoretical techniques for predicting 
aerodynamic characteristics and surface pressures for relatively slender 
vehicles at moderate hypersonic speeds were developed. Emphasis was placed 
on approaches that would be responsive to preliminary configuration design 
level of effort. Second order small disturbance and full potential theory 
was utilized to meet this objective. 

Numerical pilot codes were developed for relatively general three 
dimensional geometries to evaluate the capability of the approximate 
equations of motion considered. Results from the computations indicate good 
agreement with higher order solutions and experimental results for a variety 
of wing, body, and wing-body shapes for values of the hypersonic similarity 
parameter MS approaching one. Case computational times of a minute were 
achieved for practical aircraft arrangements. 
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1. INTRODUCTION 


An examination of the literature for airbreathing hypersonic concepts 
provides an indication of the flexibility and generality required for a 
prediction technique. Typical configuration development variables include 
wing section, incidence, height, dihedral, planform, effectiveness of 
longitudinal control surfaces for trim, effectiveness of empennage for 
directional stability, and propulsion system- airframe interactions. 

State-of-the-art response to these prediction requirements is provided 
by hypersonic impact methods as well as linearized analysis and design 
algorithms. These approaches can treat complex geometries with minimum 
response time and cost, with efficient predicted data coverage in terms 
of Mach number, angle of attack, trim deflection, yaw angle, etc. Shortcomings 
are present, however, in both the impact and linearized methods. For the 
former, interference between surface elements is totally ignored in 
implementations such as classical Newtonian, tangent wedge, and cone theories. 
Cross- flow interactions and stagnation point singularities are also implicitly 
disregarded. In the latter, shocks, vorticity, and entropy wakes and layers 
are excluded. Furthermore, superposition of elementary solutions such as 
those for thickness and angle of attack freely used in linear models are, 
strictly speaking, invalid at hypersonic speeds. 

A need exists for new aerodynamic prediction techniques to optimize 
vehicles designed to travel at supersonic/hypersonic speeds. One requirement 
of a new aerodynamic prediction is that it be more accurate than simple non- 
interfering panel methods. Another specification is that it be more computa- 
tionally efficient than currently available explicit finite-difference methods 
so that it can be incorporated into a practical design procedure. The new 
approach should include enough of the physics of the flow to allow realistic 
optimization and should permit consideration of appropriate interactions 
between components of promising hypersonic arrangements, since this has been 
found to be the key to increasing aerodynamic efficiency at supersonic speeds. 
Nonlinear potential theoretical formulations hold the promise of meeting this 
objective and providing economic design codes which are responsive to 
preliminary vehicle definition efforts. 
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2. LIST OF SYMBOLS* 


AR 

b 

c 

e 

CD 


cl 

Cm 


C P 

D 

e 

k 

L 

M 

n 

P 

q 

s 

t 


u,v,w 


a 

S 

Y 


aspect ratio, b^/S 
wing span 
chord 

mean aerodynamic chord 

drag coefficient, D 
qS 

lift coefficient, L 

qs 

pitching moment coefficient M 

qSc 


pressure coefficient, P-P°o 

q» 

drag 

unit normal 

second order coefficient, 
lift 


Mach number or pitching moment 
normal 

static pressure 
dynamic pressure 
reference area 
maximum thickness 

axial, lateral, and vertical perturbation velocity 
axial, lateral, vertical body axis cartesian coordinates 
angle of attack 


/ M^oo-1 

ratio of specific heats 
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6 flow deflection angle 

7 ] wing span station /(b/2) 

0 dihedral angle of aerodynamic surface 

X taper ratio, ct/c r 

A sweep 

<j> ' first order velocity perturbation potential 

<j) (2) second order velocity perturbation potential 

$ total velocity potential 

ip (x,y) surface displacement 


ac 

b 

LE 

r 

t 

x,y,z 

& 

u 

c,y 

t,r 

n 


a 


SUBSCRIPTS 
aerodynamic center 
body 

leading edge 

root 

tip 

derivative with respect to x,y,z respectively 

lower surface 

upper surface 

camber 

thickness 

normal 

free stream 

angle of attack 


*Additional specialized nomenclature is defined in the theoretical sections 
at the point of introduction 
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3. METHODOLOGY 


Emphasis is placed on approximate theoretical approaches which are capable 
of treating relatively general three dimensional problems but still sufficiently 
simple to be responsive to vehicle preliminary design efforts. The basic 
intent of the methodology is to produce improvement in lift-drag ratio of 
hypersonic cruise vehicles. As a result of the strong impact that favorable 
interference has had on supersonic design and the use of such concepts in recent 
advanced hypersonic aircraft studies, candidate analysis should be general 
enough to systematically treat such problems. Finally, interest in high 
aerodynamic efficiency emphasizes relatively slender configurations at modest 
angle of attack; that is, moderate values of the hypersonic similarity 
parameter. 

1 2 

Prior theoretical effort ’ has advanced the supersonic/hypersonic 
aerodynamic prediction state-of-the-art at the preliminary design level. 
Numerical second-order potential small disturbance analysis was developed as a 
first step up from the widely used linear theory. Such a formulation 
incorporates nonlinear behavior by iteration of the Prandtl-Glauert 
approximation. This approach is known to extend the prediction success for 
airfoil and cone surface pressure to substantially higher values of the 
hypersonic similarity parameter than the first-order theory. A typical 
three-dimensional prediction improvement provided by the numerical second- 
order analysis ^ is presented in figures 3.1 and 3.2. These results support 
extensions to treat multiple-surface analysis and wing-body design problems 
in the present study. The next level of theoretical richness vis-a-vis a 
full -potential equation of motion formulation was explored in parallel^. 

This analysis eliminates edge singularities and improves treatment of 
characteristic surfaces but still retains the isentropic assumption. Typical 
prediction success for a conical wing-body is presented in figure 3.3. 

These results indicate sufficient promise to pursue extension of the analysis 
to general wing-body configurations. 

Hypersonic small disturbance theory was considered in an earlier study! 
in recognition of the progressive non- isentropic behavior of the flow as the 
value of the hypersonic similarity parameter increases. Finite difference 
analysis of this approximation^ indicated that the solution was essentially 
as complex as that for the Euler equation and thus would not be particularly 
responsive to preliminary design level of effort. This approach was not 
pursued in the present study on the basis of this finding and the previously 
cited success of potential analysis at moderate hypersonic conditions. 


4 



Figure 3.1. Predicted Aerodynamic Center Location of a 
Subsonic Edge Delta Wing 


















Figure 3.3. Predicted Compression Surface Pressures for 
a Conical Wing-Body at = 4 , a = 5° 
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SECOND ORDER POTENTIAL ANALYSIS/DESIGN 


4.1 APPROACH 

The solution accurate to second order for the flow around wing-body 
combinations of aerodynamic interest may be represented by the velocity 
potential, <£ , written as the sum of the first and second order velocity 
potentials. 


<£(x,y,z) = ' 0 (x r y,z) + ^Xx,y,z) 


The condition of no flow through the boundary requires that 


[ cosa., e x + sina,,, e^ + vf ]. n * 0 


on the surface of the configuration. The solution , 0, accurate to first 
order, must satisfy the Prandtl-Glauert equation, 


. 32 32 32 

□N = i (1-Mj — + — + — r 0 = 0 
3x2 3y 2 az 2 


and the the boundary condition 


[e+a_e+ v0 3 - 7i = 0 


on the surface. For the configurations of interest the aerodynamic surfaces 
are assumed to be thin, and the first order boundary condition is satisfied on 
the mean camber line in accordance with thin wing theory. 


30 

3Y 

3n 

ax 


a w cosfl 


where, 6, is the local dihedral angle and, z = ¥(x,y), is the equation for the 
surface deflection written in a coordinate system with the z axis 
perpendicular to the local mean camber plane, and the x and y axis located 
within this plane, with the x axis in the streamwise direction. To first 
order there is no effect on the solution if a distinction is made between the 
angle of attack of the freestream and the angle of attack of the 
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configuration. To second order there also should be no difference in the 
solution, but for numerical computation purposes a distinction should be made, 
and the angle of attack should be accounted for by a B , and not by ¥(x,y). 


The method used to obtain a first order solution made use of axial 


singularities and quadrilateral source panels on the body surface in 
conjunction with quadrilateral source and vortex panels to represent the 
(thin) aerodynamic surfaces. The vortex panels, used to represent lift, were 
of constant strength, while the source panels, used to represent thickness, 
could be of constant strength or could be made to vary linearly in both the 
chordwise and spanwise directions. The vortex panel and body source panel 
strengths where obtained by solving a set of simultaneous linear equations and 
thereby satisfying the boundary conditions at a set of control points. The 
singularities used to simulate a multiple surface-body configuration are 



Distribution of Singularities on an Arbitrary Configuration 
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When the panel singularity strengths were obtained the first order flow 
properties could be determined anywhere in the field. The second order 
potential, utilizes the first order solution. It must satisfy the 
nonhomogeneous Prandtl-Glauert equation: 


■ »C~>[(l-iC)+ — M» ] f f + [0 + a„ ]* ■ 

3x 2 * a a . J 


as well as enable the boundary conditions to be satisfied to second order on 
the surface. 

A second order solution using the source volume formulation described in 
reference 1, encounters numerical difficulties for complex configurations, 
wings with subsonic edges, and supersonic wing-body configurations. The 
primary reason is that source volume strengths require the calculation of 
spatial derivatives, and these cannot be obtained accurately enough from the 
first order solution when using the panel formulation. In addition the 
velocity discontinuities present in supersonic flow are accentuated, since the 
discontinuities introduced by the panel corners and edges do not attenuate 
with distance. Also no reasonable mesh density of spatial source volumes can 
properly account for the large gradients in flow properties near subsonic 
edges. 


Therefore an approximate method which does not require the use of a 
spatial distribution of sources was developed^ This method requires no 
calculation of flow properties off of the surface, and can be used even when 
large gradients are present (e.g. with subsonic edges). The method of 
solution is a three dimensional, modification, based on an infinitely skewed 
wing transformation, of the exact second order solution for the pressure 
coefficient on thin airfoils in' two dimensions. Using this method, the 
solution for the thin aerodynamic lifting surfaces is given by the following: 


Cp(x,y) = - 2 


0 x + k(M) E^(l-M^) i + + ;{^+a J - + *4+ \ 

■* a *o 


where cf <p = 0 , k(M) = 


1 


(1-M2) 


^ 2 ■ Z 

7+1 M 4 




i 1 + 


4 ( l -M2 ) 


and <$Jx,y,0) 


(1-M J <t> ¥ + <t> Y + V [ (1-MJ 0 + 4> 3 

* * a a ** a* 
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Figure 4.1.1. Second Order Coefficient 



¥(x,y) represents the upper or lower surface of the airfoil. This solution, 
in effect, represents a local type solution, and all first order velocities 
represented above should be in a coordinate system rotated to make the z-axis 
normal to the local surface (e.g. with nonplanar configurations), but with no 
local angle of incidence. 

When applied in three dimensions the transformed two dimensional equation 
should give the exact second order pressures whenever the flow is locally two 
dimensional and may be approximated by an infinite skewed system. The Mach 
number, M, used in calculating the coefficient k(M) is equal to 

= M,, cos A (i.e. the normal Mach number) 

Whenever M = M^. is equal to unity, k(M) becomes infinite, indicating the 
second order solution is not valid near M * 1. If M. is used for calculating 
k(M) the solution remains finite except near » 1, Since k(M) is a slowly 
varying function of M (see page 11), except near M» = 1, Mm,,) can be 
substituted for k(M^) without a significant change in the result^ 

The second order modifications to the pressure on the body are 

influenced only by the function used to satisfy the second order boundary 
condition. However, in a sense, the second order requirements on the body are 
not as stringent as on the aerodynamic surfaces, since the boundary conditions 
on the body are satisfied on the actual surface, rather than the mean camber 
line. In addition the nonlinear expression for Cp in terms of the velocities 
is used on the body. The terms missing from the complete second order 
solution on the body are due to the spatial term in the non-homogeneous 
Prandtl-Glauert equation for 


4.2 FULL CONFIGURATION ANALYSIS 


The first order solution is obtained by solving for the panel strengths 
required to satisfy the first order boundary conditions. The panel strengths 
enable the calculation of first order perturbation velocities, used for the 
calculation of second order velocities, at panel control points. Since the 
first order solution is linear, the solution and the resulting perturbation 
velocities can be considered to be the appropriate sum of the following four 
basic solutions: 

1. Contributions from all panels due to an incremental 

change in angle of attack (add load velocities). 

2. The incremental contribution from all panels due to 

the twist and camber of the aerodynamic surfaces. 

3. The velocity perturbations due to thickness. 

4. The velocity perturbations due to the body in axial flow. 
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On the thin lifting surfaces, (chord plane vortex panels), some 
components have an odd or even part with respect to the local z ■ 0 plane. 

Odd symmetry terms are those which are of equal magnitide but of opposite sign 
across the z ■ 0 plane, while even symmetry terms are those which are 
continuous across the z ■ 0 plane. Therefore if we define: 


a ■ angle of attack of the free stream. 

7 ■ a camber coefficient ( ■ 1.0 for camber, * 0.0 without ) 

r « thickness coefficient ( * 1.0 with. = 0.0 without) I 

b = body contribution 


Then we can write, 

u a $ » o ( u ± ?Au ) 

v » $ » a ( v ± jAv ) 

w * 6 » aw 


+ 7 ( u y ± jAu t ) 
+ 7 ( v^± jAv,) 
+ n, 


+ T ( U T ± £Au t ) + ( u t ± £Au t ) 

+ T ( V T ± |Av t ) + ( v t ± jAvJ 

+ r ( ± \) 

* 


where ± indicate upper or lower surface. The contributions having a ± have an 
odd symmetry with respect to the local z = 0 plane. 

The surface displacements V are composed of three parts. 

1. Displacements due to angle of attack of the configuration. 

2. Displacements due to twist and camber. 

3. Displacements due to thickness. 

Therefore we can write, 


¥ *» a ¥ + 7 ¥ i rf 

« * T 

All contributions to the second order terms are composed of products of first 
order terms. These terms will have either an odd or an even symmetry with 
respect to the local z « 0 plane. The odd symmetry terms in the Cp expression 
contribute to a ACp across the aerodynamic surface while the even symmetry 
terms contribute to an equal Cp on each side. Therefore since, 

Cp = - 2 k(M)| + iU/ «J -|«i + <p, ■ 
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where the boundary condition for the potential function c$>(x,y,z), used 
to satisfy the second order boundary condition is, 

? (x,y,0) = (1-0* * +♦*+*[ (1-0* + * 1 

with 

tt*$(x,y,z) * 0 


we can write the expression for the second order contribution to ACp 
as: 


<*.) f 

ACp = - 2 k(M)j(l-M2) [ou^ + 7 U X + TU r + uJIaAu^ + 7 Au^ + rAU T + AuJ 
+ [av„ + 7 v„ + rv + v ][aAv + 7 Av + rAv + Av ] 

• Z b M Jf lb 

♦ 2 77 [“A/ * VrJ + A \ } 


where 


7U x + 

TU + 
T 

u b > + 

ir%.( aAu + 7 Au + rAu + 
2 r * ' T 

AuJ 

7v^ + 

7V_ + 
Z. 

v K ) + 

±T% (aAv + 7 Av + rAv + 

h n * r 

i» k )] 

7U y + 
* * 

ru, + 

u w ) + 

( av + tv + tv + 

X 

v» 


+ ir¥J(l-M*)(aAu, + 7Au„ + rAu + Au. )+ (aAv + 7Av + tAv_ + Av )] 

4 T "x ** \ \ 3 A ‘3 S 

The expression <$ x is calculated by differentiating the solution for the 
function which is used to satisfy the boundary conditions to second order. 
The odd symmetry second order boundary conditions will be satisfied by a 
solution using source (thickness) panels, while the even symmetry second order 
boundary conditions will be satisfied using vortex panels and will contribute 
to a net ACp across the surface. 

By setting the various coefficients a, 7, or r equal to zero, the second 
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order effect of various combinations of angle of attack, camber ( 7 ■ 1.0 ) or 
no camber ( 7 » 0.0 ), and thickness ( r * 1,0) or no thickness ( r « 0.0 ), 
may be analyzed from the basic solutions. The term Au*. comes from the lift 
generated by the thickness and is generally small, and in fact is zero for 
planar configurations. The term Au b constitutes the lift induced by the body 
in axial flow and is also, in general, small. For planar configurations or 
when there is symmetry about the z * 0 plane, the. velocity expressions 
simplify, 


u = 


i ■ a ( u ± rAu ) + 7 ( u + rAu ) + r u- 

< at** Q * * Z 


* a ( v^ ± £av ) + 7 ( v v ± jAv y ) + r v 


a 


+ 7 w. 


± jrA'F r 


4.3 WING-BODY DRAG OPTIMIZATION 
4.3.1 Force Calculation 


The first or second order formulations for the near field force 
determination for any configuration are similar. Basically the calculation 
involves multiplying the local value of Cp times the projected area. The lift 
and drag coefficients can be written as integrals of the pressure coefficient 
over the surface area. 


C L 


C D 


Cp n z dS 


J Cp n x dS 


where the local normal is n = (n ,n ,n ) 

X il * 


In general for thin aerodynamic lifting surfaces, 


drag/area * ( Cp^v^ - Cp^ ) 

= - 0.5 ( Cp - Cp ) ( w + w ) 

X u. Vfc 

+ 0.5 ( Cp + Cp )(»“*.) 
= - ACp w c + 2 Cp w 

where u and X refer to upper and lower surfaces, and: 


wv = 

0.5 ( V - \ ) 

(camber normal velocity) 

wt = 

0.5 ( w^- w^ ) 

(thickness normal velocity) 

ACp = 

( Cp - Cp ) 
X 

(normal force Cp) 

2 Cp - 

( Cp + Cp ) 

X U 

(even symmetry Cp) 
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Using this formulation the integrals may be changed to a sum over all 
vortex and thickness panels, which are used on the aerodynamic surfaces, and 
over all source panels, which are used on the body. 


C D 


- 2 


NTV 




NTB 



i=l 

ACp. 

4 

COS0. 

4 

Pa < ' 

- 

i=l 

Cpb. nz^ 

Pab. 

4 

NTV 




NTB 




ACp. 

WV. 

Pa. 

- S 

Cpb. nx. 

Pab . 

i-1 

i 

4 

4 

i»l 

4 4 

4 

NTV 








cp. 

wt. 

Pa. 




i=l 

4 

4 

4 





where 


ACp. * 

4 

the 

Cp. ■ 

L 

the 

Cpb : = 

the 

WV. » 

L 

the 

wt ; » 

the 

COS0. a 

L 

the 

nx . = 

4 

the 

nz . = 

i 

the 

Pa . a 

4 

the 

Pab. = 

the 


On each panel the odd symmetry normal velocity, wt, is due to the 
thickness, while the even symmetry normal velocity, wv, is primarily due to 
camber (or angle of attack). The total number of vortex panels used to 
represent the configuration is equal to NTV, and the total number of body 
panels is equal to NTB. The Cp's on the vortex (or thickness) panels may be 
further divided into Cp's of order one and two. 
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ACp 

Cp 


o» u> 

* ACp + ACp 


to «) 

3 Cp + Cp 


The Cp on the body, Cpb- t , can be written in terms of the induced 
velocities. To first order: 


Cpb L = - 2 Dx . 


where the perturbation velocity at the control point of panel i is: 


U. - (Ux.,Uy. ,Uz.) 

«• i u <■ 


and 


U.» n 


Un. 

L 


If desired, higher order formulas may be used, although in the present 
formulation the velocities induced on the body are not fully accurate to 
second order. 

In performing the optimization we can assume the thickness distribution 
and body geometry to be given. The optimization procedure consists of finding 
the body angle of attack and, the twist and camber of all lifting surfaces. 
This optimum results in a minimum value for the drag expression, while 
satisfying a set of constraints on lift coefficient, center of pressure, etc. 


4.3.2 First Order Optimization 


To first order the thickness drag is nearly independent of lift, and 
therefore may be neglected in any optimization procedure. Therefore only the 
induced ACp's on the vortex (thickness) panels must be considered. These 
panel ACp's will be simply written as Cp^ . For a fixed body geometry, the 
local angle of attack at the vortex panel control points, Un^, and the body 
angle of attack « , may comprise a set of independent variables. Then a set 
of dependent variables which may be solved for are: 


1. The vortex panel strengths Cp . j=l,NTV 


2. The source panel strengths < j . j=l,NTB 
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where NTV and NTB are the total numbers of vortex and source panels 
respectively. The equations used for the solution are the conditions of flow 
tangency at the vortex and source panel control points. Using the computed 
aerodynamic influence coefficients, A, for each panel we can express the normal 
velocity on each panel in the form: 


NTV NTB 

Un. = T An..Cp. + Anb. a i = l,NTP 

L pL J A-l ^ A 


where NTP = NTV + NTB 


Alternately, the source strengths cr, can be solved for as a function of the 
Cp's on the vortex panels, and a. 


NTB 

NT V 



<ET Anb. a = 


An.. Cp - nz. a 

i=l ,NT3 

51=1 tSL * 

j-1 

<■4 1 



NTV 



°x 

j-1 

B Cp - B a 

j 

£=1,NTB 


where m = NTV + 1 

This means we can consider the independent variables to be the Cp's on 
the vortex panels and a: 

(Cp., a ) j = 1, NTV 


Everything of interest may be written in terms of these variables. For 
example, using the aerodynamic influence coefficients for the x-velocity, we 
can express the x-velocities in terms of the dependent variables: 


Ux . 


L 


becomes 



Ax..Cp. + 
j 


NTB 

"21 Axb. cr 
o=i ‘A 


i = 1 , NTP 
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NTV 


NTB 


NTV 


Ux . 

s 

21 

Ax.. Cp + 

Axb. 

21 

B . 

L 



‘4 J 


c4 | 

_ j-i 

-S-4 



NTV 

NTB 

- 



NTB 


s 

z: 

Ax ..+ ]>“ 

L LJ 

Axb. B 
*•4 Aj 

Cp 

4 

i 

z: 

s-l 



NTV 







= 


Bx . . Cp + 

Bx. a 






j 3 l 

*•4 j 





since 

on 

the body panels, 

Cpb^ = 

- 

2 Ox 

i 


iJL 


NTV 

Cpb. * > Bp,.Cp. + Bp a 
j=*l *-i J i~*. 

and 


NTV 

NTB 

3 > Cp. cosfi. Pa. - 

> Cpb . nz Pab . 

i 3 1 *• L 1 

i 3 1 *■ 4 4 


becomes 


NTV 


= > 

Cp cosfl. Pa. 

i=l 

i U L 

NTB 

NT V 

- ">~ 

nz. Pab. ( > Bp Cp + Bp a ) 

i=l 

fc *• 4 L~. 
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The drag coefficient may be written as a quadratic function of the set of 
independent variables (Cp. , a ). 


NTV 



NTB 



becomes 

NTV 

Cn = - >• Pa . Cp. 
i-1 L c 


NTV 

>• 3n..Cp. + Bn. a 

“fZT ‘--I -i 



Q a 
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where n = m + 1 = NTV + 2 , and 



Q. = 


Q. = 

L'W 

Q = 


Pa. Bn . . 

<. *-a 

Pa. Bn • 

I «•— ■ 

NTB 

21 n Pab Bp 
k=i * a . -*■ *<• 

NTB 

21 n Pab Bp 
k=l *k * 


Other constraints may be written as a linear function of the independent 
variables (Cp- , a ) in a similar manner. The optimization consists of 
minimizing the expression for Cd subject to the constraint equations. Since 
the expression for Cd is quadratic in the unknowns, the minimization may be 
performed using the method of Lagrange multipliers. The unknowns can then be 
used to solve for the Cp's and also the normal velocities on the aerodynamic 
surfaces. The normal velocities on the aerodynamic surfaces determine the 
twist and camber. 
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4.3.3 An Exact Second Order Optimization In Two Dimensions 


In two-dimensional supersonic flow an exact 2nd order solution is 
available. This example will serve to demonstrate the necessity of including 
the thickness drag in optimizing drag due to lift, when second order pressures 
are considered. 


* 

— ( w + w ) + 2 k 

/a * «■ 

cp * “ 

— ( w - w ) + 2 k 

p * «■ 


2 . 2 . 

where . J2> - M,,, - 1 

and k ** MM.) where, k(M) 


(1-M2) 


1 + 


7+1 


M4 


(1-M2) 


therefore 

ACp 
2 Cp 

and drag/area 


^ 1 
- [ w + 2 0 k w w ] 

pc. 


c : t 


- - [ w + 0 k ( w +w z )] 

p * t c 

= - [ w + 2 0 k w w ] w„ 

p> c <it c 

+ - [ w + J k ( w*t w l ) ] w 

p> t t c. i 


= ~ [ ( 1 + 3 0 k w ) w 2- + ( 1 + 0 k w ) w‘] 

/ s> t c t t 

Therefore the thickness and camber interact to give additional drag terms from 

both 


- ACp w e and 


2 Cp w^ 
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With a fixed thickness distribution the drag may be minimized subject to 
a given lift. Ignoring the terms which depend on thickness only: 


minimize 


CD 




I 

( 1 + 3 0 k w ± ) w c * dx 

« 


subject to C 


^ ■ ■ ji 1 


(1+ 2 0 k w ) v dx 


The solution may be obtained using a variational approach. 


let * c (x) = W(x,X) + 5fl(x) 


where 5n(x) is an arbitrary function which is allowed to vary, and X 
is a constant which can be adjusted to satisfy the C L constraint. If 
W(x,X) is the function which gives the minimum drag, subject to the C L 
constraint, then the addition of any function 6n(x) will increase the 
drag. Therefore: 


~ (C D - X C L ) = 0 at 5 = 0 


and since — w,(x) = ^(x) 

c 

J [ 2 { 1 + 3 0 k w^ ) W X{ 1 + 2 0 k w c ) ] U(x) dx = 0 


Since this relation holds for all functions T(x), we must have: 


W(x,X) = 


X [ 1 + 20k w t { x ) ] 

2 [ 1 + 30k wjx) ] 


but since 


4 X 
0 2 


[ 1 + 20k w t ({) ] 
[ 1 + 30k w t ({) ] 


d* 


4 

0 


f [ 1 + 20k w t (5) ]' 
[ 1 + 30k w t ({) ] 


d{ 
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or 


2 



[ 1 + 20k v*(x) ] 
[ 1 + 30k v fe (x) ] 


0 

then w (x) « - - C L 

c 4 ■ 

[ 1 + 20k v.({) ] * 



J [ 1 + 30k v t (£) ] 


and therefore 


C D 


OPT 


0 z 

~ C L 

4 


[ 1 + 20k w t (J) ] 
[ 1 + 30k w t (£) ] 


d{ 


The drag value is slightly lower than the result using a first order 
analysis. The first order result is obtained by setting k * 0. 


Cjjoor 

c* 


0 

4 


The second order drag may be evaluated for any specific thickness 
distribution. As an example, a biconvex thickness distribution can be used. 


z t (x) = 2 r x { 1 - x ) vr t (x) =2 r ( 1 - 2 x ) 

First consider the following: 


{ 1 + 2 a ) 
{ 1 + 3 a ) 


*■ , *• 
2 ( | + 3 a ) 

3* ( 1 + 3 a ) 


1 [ j + ( 1 + 3 a ) ] ^ 

3* ( 1 + 3 a ) 


z 

2 



1 

• + 

2 * { 1 + 3 a ) 


{ 1 + 3 a ) \ 
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Now let 


a(x) = 0 k wjx) = 2 0kT ( 1 - 2 x ) 


= — e ( 1 - 2 x ) 
3 


where e * 6 0 k r 


Then let 
f 


dx 


( 1 + 3 a ) 


dx 


[ l + e ( 1 - 2 x ) ] 


1 1 

- - log [ 1+e ( l-2x ) ] 
e 2 


1 1 (1+e) 

- - log 

e 2 (1-e) 


Therefore: 


[ 1 + 20 k W t (x) ] 
[ 1 + 30 k w 4 (x) ] 


dx 


1 

9 

1 

9 


O 

J 


( 1 + 3 a ) 


( 1 + 3 a ) 


[1+2 a(x) ] 
[1+3 a(x) ] 


dx 


+ 8 + 12 a 


dx 


+ 8 + 4 e ( 1 - 2 x ) 


dx 


1 

( f + 8 ) 

9 
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0 

and C n = - 

u ° p T 4 

Since lim f = 1 

e -» o 

order result as e -»■ 
If we set w c (x) = c 


C L 

( f + 8 ) 

this result reduces to the first 
0 

(the first order result), we get 


C D 

C L 



+ 3 0 k w^(x) ] c Z dx 

c 

+ 2 0 k w t (x) ] c dx 


O 

and whenever the thickness distribution integrates to zero. 


r 

Dplat 



J [ 1 + 3 0 k w ± (x) ] dx 

A O 



+ 20k w fc (x) 




The ratio of the optimum camber distribution drag to the uncambered 
airfoil drag at the same value of C L is: 


CDoi» r 

The angle of attack of the optimum camber airfoil also differs from 
first order result, since there is some lift due to the camber. The 
angle of attack is the integral of w t (x) over the chord. 

i 

a = [ w c (x) dx 


9 

( f + 8 ) 
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therefore 


CV 

OPT 


[l+-|e(l-2x)] 

2 dx 

J[l+ e! l*2x)] 

0 

- c L 

4 j ( f + 8 ) 


( f + 2 ) 0 

3 - C L 

( f + 8 ) 4 


The results for the biconvex airfoil will be presented and discussed in 
section 6.1.2. 


4.3.4 Second Order Wing Optimization 


Again the optimization problem consists of mimimizing the expression for 
drag while satisfying a set of constraints on C L , center of pressure, etc, and 
again we will assume the thickness distribution to be given. However if 
second order pressures are considered, the influence of the lift on the 
thickness drag must be considered. This was demonstrated in the previous two 
dimensional example. Since in the calculation of second order Cp's, 
derivatives of camber geometry must be computed, the camber optimization is 
conducted by considering the camber to be composed of a sum of unit camber 
functions with smooth derivatives. The set of functions must be complete 
enough to describe any optimum. There would be a set of functions at each 
span station consisting of a local angle of attack (twist) and some additional 
functions which integrate to zero in the chordwise direction. 


Let zv(x,j,k) be the camber z/c value, at x/c = x, for the k'th camber 
function, on span station j. An appropriate set of camber constraint 
functions might consist of the following functions for each span station j. 


z/c 


x 


j wv(x) dx 

o 


0 for all k, at all span stations # j 


i.e. the camber constraint function at span station j contributes 
no camber at any other span station. On span station j: 
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z/c a 

zv(x, j , 1 ) 

K — 

X 




= 

zv(x, j ,2) 

S 

x(l-x) 




a 

zv(x, j ,3) 

B 

x(l-x) 

(1- 

*2x) 


s 

zv(x, j ,4) 

a 

x(l-x) 

(1- 

-4x) 

(3-4x) 

s 

zv (x , j , 5) 

s 

x(l-x) 

(1- 

-2x) 

(x-A) (x-B) 

where 







A = 

( 2 -/2 ) 

/ 4 





B = 

( 2 +/2.) 

/ 4 





Therefore 

the camber 

and i 

: amber : 

normal 

velocities for 


written as the matrix product of a constraint matrix (Zjj* for the 
displacements, zv, and Cjj^for the normal velocities, wv), and a vector 
of the unknown coefficients, x A , of each camber constraint. If N is the 
total number of constraint functions: 


zv. 

t. 


N 

k-1 


Z x 


i = 1, NTV 


WV. a 

L 


N 

k=l 


C 

tk At 


i » 1, NTV 


A set of unit solutions for the first order Cp's result from the 
camber constraint functions. This means the first order Cp's may be 
written as a linear sum in terms of the unknowns x^. If P^represents 
the first order ACp at panel i due to camber constraint function k, 
then we can write: 


Si CO 

ACp. * T P, x i - 1 r NTV 


The second order Cp's may be written in terms of the first order 
solution. 
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If we let Z(x,y) = Zv + zt 


= f'U,?) - ^(x-x ) 

+ + J(^+ a j - + Z 0^ + <? ) 

with the second order boundary condition on 

= (1-Mj 0 Z + 0 Z + Z [ (1-Mj 0 + 0 ] 

2 x * a a ** 

where Q 1 ef> = 0 


( 3 .) 


Cp = - 2 k(M)[^(l-M^)02 


The second order Cp’s may be broken into components with odd and even 
symmetry about the z = 0 plane of each panel. 


0^ * u( vortex-even) + u(thickness) ± u(vortex-odd) 

0 = v( vortex-even) + v(thickness) ± v(vortex-odd) 

* 

0 = w(vortex-even) + w(thickness-even) ± w(thickness) 

2 


The planar case is easier and will be considered first. In this case 
u(vortex-even) = v(vortex-even) = w(thickness-even) = 0 


and 

n 

X 

o 

ut 

+ 

uv 


0, - 

vt 

+ 

vv 


K ■ 

wv 

± 

wt 


2 


where the v and t components are entirely due to vortex and thickness 
panel perturbations respectively. The additional terms are usually 
small and can be added when the nonplanar case is considered. 


Planar Case 

For the case of configurations confined to a single plane the first 
order thickness and vortex solutions are decoupled. The solutions may 
be performed independently and the perturbation velocities have either 
an odd or even symmetry. It will be shown that the planar case allows a 
direct solution for the camber coefficients x v even when second order 
Cp's are used. The nonplanar case will require an iterative solution. 
Written in terms of the first order vortex and thickness velocities at 
each panel: 
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«) 

ACp = -2 kC’i) 


£ 

0 UV Ut + VV vt + (WV + <!„,) wt 

•oo 

wv Zt + wt Zv + d> ] 

* x T * 


(Z) 

Cp = -2 k(M) [j/S^uv + 
+ i(wv - 
+ wv Zv 

Where <9 satisfies □*<£ = 0 , and 
is given by: 


ut 1 ) + £(vv z + vt 1 ) 

o + K 

+ wt Zt + <P * V< ' 4 ] 

x * 

the second order boundary condition on 


<9 


OOD 

Z 


p ( uv wt + ut wv ) + 
+ Zt (0* — UV + — vv 


vv ^ Zt + vtr Zv 
3a 

) + Zv (jJ 2- — ut + — vt ) 


<9 


tMflO 

z 



uv wv + ut wt ) + 

ZV T 7 UV + — vv 
a* 3^ 


vv — Zv + vt — Zt 
aa as 

) + Zt (0* — ut + — vt 

a* ^ 


) 


Using a coordinate transformation at each panel: 


£ = x-Ty where T is the tangent of the panel sweep angle 
V = y 


we can write y derivatives in term of x derivatives and derivatives 
along lines of constant percent chord, ^ ( $ =» const) 


3 3 

3X 3f 
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Therefore we can write: 
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The vortex solution velocities can be expressed in terms of the unit 
camber solutions. Since the thickness solution is given, the second order 
ACp's can be expressed as a linear combination of the unknown camber 
coefficients x, . 

M. 
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The second order even Cp solution can be expressed as a quadratic 
function of the unknowns x . 
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Therefore the drag coefficient may be written: 
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By using the method of Lagrange multipliers, this quadratic form may be 
minimized in conjunction with the constraint functions which can be expressed 
as a linear combination of the unknown camber coefficients x, . 

M. 
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5. FULL POTENTIAL ANALYSIS 


The full potential equation either in conservative form or 
nonconservative form is frequently used for the treatment of transonic flow 
fields, where the local Mach number in general does not exceed ~1.4. 

However, if the assumptions of irrotationality and isentropicity are 
reasonably valid, then the full potential equation is expected to yield 
results comparable to Euler equations, even for supersonic/hypersonic flow 
fields. For preliminiary design studies, where quick turnaround is desired, 
the full potential methods can be an attractive substitute for expensive 
Euler methods and less accurate linear theory methods. 

A nonlinear aerodynamic prediction technique based on the full potential 
equation in conservation form has been developed, during this contractual 
effort, for the treatment of supersonic flows. A detailed description of the 
method is already presented in two published papers, which are enclosed in 
the Appendix section for convenience. Reference 4 describes the method from 
the concept of positive artificial viscosity and diagonal dominance. 

Reference 5 presents the method in a more elegant way, in terms of the theory 
of characteristic signal propagation. Since the Appendix section clearly 
describes the full potential method, only the add-on material (not included 
in the published articles) will be presented in this section. 

5.1 GRID GENERATION 

Special procedures are necessary to generate a grid system around sharp 
leading edges of wings. A typical cross-sectional view of such a sharp 
leading edged wing-body is shown in Figure 5.1. If one generates a grid 
system without any special treatment, most likely the grid solver routine will 
fail to converge near the sharp leading edge. Even if convergence is 
achieved, the resulting grid will not have a smooth distribution due to 
collapsing grid cells near the leading edge. To avoid this, the code in its 
present form uses the following procedure. Referring to Figure 5.1, a user 
prescribed boundary (line A-B) is chosen such that it divides the whole domain 
of computation into a lower domain and an upper domain. Along the curve A-B, 
the user also prescribes the desired grid point distribution. Then, the 
grid generation routine is used to separately generate the grid in the 
upper and lower domain. By maintaining a reasonable mesh near the leading 
edge, the user can avoid nearly vanishing Jacobians arising due to collapsing 
grid cells. Figure 5.2 shows an example grid generated for a sharp leading 
edge wing-body configuration. 

Another modification to the grid generation part of the co^e is the option 
of prescribing the geometry in the crossflow plane (constant r, plane) at 
discrete points as shown in Figure 5.3. Then a cubic spline routine is used : 
to distribute the grid points along the discretely prescribed geometry surface. 
Figure 5.4 shows the grid pattern generated around a typical wing-body cross- 
section whose shape was discretely prescribed. Use of this option in the 
code allows the user to input any arbitrarily shaped wing-body geometry. 
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Figure 5.1. Double Domain Grid Generation Procedure 
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Figure 5 . 2 . Grid Generation for a Sharp Leading Edge 
Wing-Body Configuration 


□ PRESCRIBED GEOMETRY POINT 

• GRID POINT LOCATION OBTAINED 
FROM CUBIC SPLINE ROUTINE 



Figure 5.3. Discrete Geometry Input Option 
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Figure 5.4. Grid Generation From Discrete Geometry Input Option 
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Another attractive option in the code is the ability to go from a 
crude (m x n) grid to a fine (p x q) grid calculation using appropriate 
interpolation routines that transfer the solution from (m x n) to (p x q) 
grid system. This capability is essential to handle wing-body combinations 
of the type shown in Figure 5.5. The calculation will proceed from the nose 
of the body to the plane where the wing originates from the fuselage using 
a (m x n) cruide grid distribution. Before starting the calculation over 
the wing the solution from this grid will be interpolated onto a (p x q) 
grid with clustering of points near the wing surface as shown in Figure 5.5. 
In this way, the shock that originates at the wing leading edge can be 
adequately resolved. 

5.2 CROSSFLOW SWITCH 

The code in its present form implements a robust density biasing 
scheme that is more continuously activated across the crossflow sonic line, 
than the previous method reported in Reference 5 (see Appendix section) . 

The density biasing as reported in equations (15) and (16) of the above 
mentioned AIAA Paper, uses a multiplier v = y (1 - a 2 /q 2 ), with y = 0 for 
elliptic crossflow and y = 1 for hyperbolic crossflow. Thus, across the 
crossflow sonic line, the quantity v has a discontinuous behavior since the 
total velocity q is always assumed to be greater than the speed of sound, a. 
This discontinuous behavior of the density biasing creates spurious 
oscillations near the shock wave as shown in the pressure contour plot of 
Figure 5.6 for cone at angle of attack. To avoid this oscillation, the 
current version of the density biasing uses the following: 


v = y 1 - 


22 c 

V 2 


with 


y = 0 when a. 


22 


y2 

a 2 


> 0 


= 1 when 3 l ^2 r < 0 


y2 

2 


V 2 

Since the quantity ^2 * ~ g° es to zero at the crossflow sonic line, the 

current version of the density biasing has a continuous behavior when it is 
turned on. The implementation of the new density biasing scheme produces a 
monotonic behavior across the shock wave, as shown by the pressure contour 
plot of Figure 5.7. 
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Figure 5.5. Proper Grid Arrangement for Wing-Body 
Combinations Using Interpolation Routines 
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6. RESULTS 


6.1 SECOND ORDER 

The wing-body small disturbance analysis2 was extended in section 4.Z 
to consider multiple non planar lifting surface problems in order to model 
more realistic aircraft arrangements. The impact of a canard on the wing- 
body configuration of figure 3.2 will be presented in section 6.1.1 to 
illustrate this capability. 

Nonlinear small disturbance optimization for a specified thickness 
was developed in section 4.3 to support advanced aerodynamic design 
studies. The analysis is sufficiently general to permit constraints on 
lift, trim, twist, and camber to be imposed. Minimum drag due to lift 
solutions for airfoil, wing, and wing -body model problems will be presented 
in sections 6.1.2 and 6.1.3. These results are the first supersonic non- 
linear potential optimum published to the knowledge of the contractor. 

IBM 370 case time for the second order analysis of 30 CPU seconds for 
a typical aircraft arrangement and paneling density indicates the solution 
is responsive to preliminary design level of effort. 

6.1.1 Multiple Surface Analysis 

Prediction of the effect of a trapezoidal canard on the delta wing-body 
configuration of figure 3.2 was evaluated using first and second order 
analysis for the finite element model of figure 6.1.1. Wing chordwise net pres- 
sure distributions due to angle of attack at M = 6 are presented on figure 6.1.2 
at four inboard span stations strongly influenced by the canard. The second 
order prediction for both the canard on and canard off cases exhibit a more 
forward chord loading which is consistent with measurements (e.g. figure 3.1) 
indicating a forward shift of the aerodynamic center location as the Mach 
number increases into the hypersonic range. The principal impact of the 
canard is to decrease the wing sectional loading inboard and increase it 
outboard of its tip location n = .25. 

6.1.2 Two Dimensional Optimization 

Minimum drag airfoil results were developed as a precursor effort to the 
three dimensional problem. An exact second order solution exists for this 
class of problems and is given by 
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The resulting optima (see section 4.3.3) for specified thickness are 
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These equations may be evaluated in closed form for an analytical thickness 
distribution. The solution for a biconvex airfoil = Ztx(i-x) is 


<*. a p r _ 3 ( ~f 2) 

Cl i- ( -f + 8) 


a Z c 
a x 


^ c - 2 - 

+ " (i+8) 


[ I + f<=Cl-2.x>] 

[ 1 + <e. c » - z>0] 


_ _£•_ «? 

C * ~ i- (•¥■+ 8) 


where 


e 


< I 


44 • 


4 


i -e 



The prediction for this case is presented on figure 6.1.3 as a function of 
the second order similarity parameter 8 k t which accounts for the nonlinear 
coupling between thickness and lift. First order (i.e. linear) levels are 
recovered in the limit e -* 0. Examination of the results indicate the 
optimum nonlinear drag levels become progressively lower than standard 
linear results as the Mach number or thickness or both increase. The 
reduction requires a small positive camber as indicated on the bottom of the 
figure. Nonlinear small disturbance compressibility also produces a 
more forward aerodynamic center location which is consistent with the three 
dimensional results of section 6.1.1 and figure 3.1. 

6.1.3 Three Dimensional Optimization 

First and second order trimmed optimum drag due to lift for the model 
problem of figure 6.1.4 is compared on figure 6.1.5 as a function of the 
longitudinal stability parameter dCM/dCL at M = 2 and a design Cl = 0.1. The 
abscissa scale may be alternately interpreted as the pitching moment at 
zero lift % = -CLp dCM/dCL. The second order nonlinear solution considers 
the coupling between thickness and lift for a five percent biconvex airfoil. 
Uncambered (i.e. flat plate) linear and nonlinear results are also presented 
to judge the potential benefit of optimization. Examination of figure 6.1.5 
indicates a five percent improvement in lifting efficiency is achieved 
relative to the best linear optimum and occurs at approximately twice the 
longitudinal stability level (dQn/dCL = -.072 versus -.14). The three 
dimensional supersonic nonlinear optimum results are the first published 
to the knowledge of the contractor. 

The first and second order optimum twist -and camber for the minimum drag 
points of figure 6.1.5 are presented on figures 6.1.6 and 6.1.7 respectively. 
Both solutions have been constrained by the use of a linear combination of 
analytic functions to produce smooth deformation. 

Second order analysis of the first order optimum and first order analysis 
of the second order optimum are also presented on figure 6.1.5 to establish 
the relative importance of the pressure and deformation differences between 
the linear and nonlinear solution. The former consideration is the more 
important of the two for the case under consideration. The results also 
provide insight into the test finding that linear designs typically exhibit 
maximum drag due to lift efficiency at a lift coefficient higher than the 
design condition. This is apparently due in part to the reduction in drag 
due to lift factor (for stable balance) resulting from nonlinear 
compressibility effects. 

Optimization was also performed for a supersonic edge condition at 
M = 4.0. The results are presented as figure 6.1.8. The impact of twist 
and camber is small for both the first and second order computations. The 
former result is well known and has led to the use of planforms having a 
subsonic edge at cruise. 
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Figure 6.1.1. Wing-Body-Canard Finite Element Model 
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Figure 6.1.2. Effect of Canard on Wing Chord Load Due to Angle of Attack at M 



















Figure 6.1.3. Second Order Optimum for a Biconvex Airfoil 
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Figure 6.1.4. Three Dimensional Second Order Optimum 
Model Problem 
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Figure 6.1.5. First and Second Order Optima for a 
Subsonic Edge Condition 
M = 2.0, C L = 0.1 
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Figure 6.1.6. First Order Optimum Solution for dCM/dC L 































Figure 6.1.7. Completed 
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Figure 6.1.8. First and Second Order Optima for 
a Supersonic Edge Condition 
M = 4.0, Cl = 0.1 


55 











6.2 FULL POTENTIAL 

Results from the conservative full potential formulation are presented 
for both conical and nonconical supersonic flows. The published papers 
enclosed in the Appendix section cover results for a wide variety of 
configurations that include a supersonic delta wing, cones at angle of attack, 
conical thin elliptic wings, and wing-body combinations, circular arc- 
cylinder body, and nonconical arrow-wing body combination. 

One additional result not presented in the Appendix is reported in 
this section. Figures 6.2.1 and 6.2.2 show the grid arrangement at 
various inarching planes for the Sears-Haack wing-body combination. The 
grid is smooth around the sharp leading edge generated using the double domain 
grid generation procedure described in Figure 5.1. The pressure distribution 
on the flat upper surface of the wing for = 6, a = -4° is shown in Figure 
6.2.3 at the axial station x/l = 0.759. These are preliminary results from 
the code and further improvements are possible with the use of the newly 
developed discrete geometry input option. This will be pursued in the 
follow-on program. 
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Figure 6.2.3. Sears - H aac k Wing-Body Pressure Distribution at 
Moo 3 6.0 and a =-4° 
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7. CONCLUSIONS 


Based on the theoretical development and comparison with higher order 
results /experimental measurements described in this document, the 
following conclusions are made. 

1. Improved prediction of supersonic/hypersonic aerodynamic 
characteristics and surface pressures for general wing-body 
shapes has been demonstrated using nonlinear potential 
analysis for values of the hypersonic similarity parameter, 

M6, less than one. 

2. Second order theory provides a systematic means of extending linear 
analysis for full configurations. Thirty second CPU solution 

time /Mach number is typical for a wing-body-canard analysis or 
wing-body optimization problems. 

3. Full potential analysis successfully eliminates subsonic edge 
singularities and linear characteristics approximations of 
second order theory. The formulation is an order of magnitude 
faster than Euler solvers while maintaining comparable 
prediction accuracy. 

4. Potential theory provides an advanced aerodynamic prediction 
technique that is responsive to the preliminary design problem 
at moderate hypersonic conditions. 
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for Supersonic Flows 
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An aerodynamic prediction technique based on the full potential equation in conservation form is developed 
for the treatment of supersonic flows. This technique bridges the gap between simplistic linear theory methods 
and complex Euler solvers. A local density linearization concept and a second-order-accurate retarded density 
scheme, both producing the correct artificial viscosity, are introduced in developing an implicit marching scheme 
for solving the scalar potential Results for conical flows over delta wings, cones, and wing-body com- 
binations, and for nonconicai flows over bodies of revolution at angles of attack are compared with Euler and 
nonconservative full potential calculations and experimental data. The present formulation requires an order of 
magnitude less computer time and significantly less computer memory. over Euler methods. 


I. Introduction 

A ERODYNAMIC prediction techniques that can handle 
significant geometric complexity for use in supersonic or 
hypersonic configuration design are based on either hyper- 
sonic impact methods' or linear theory analysis, 2 both of 
which require minimum response time and cost. However, 
shortcomings are present in both the impact and linearized 
methods. Aside from these simplified techniques, limited 
capabilities also exist for calculating supersonic fiowfields 
using very complex Euler codes, 3-4 using either shock cap- 
turing 3 or shock fitting* 4 methods. The use of these codes as 
viable aerodynamic prediction techniques for configuration 
design is, however, not practical due to their slow response 
time (requirement of large computer memory) and excessive 
computer cost per run due to strict stability requirements. 
Thus, we have on one end of the spectrum, very simplified 
codes that require minimum computer time to provide less 
accurate results and, on the other end, very complex Euler 
codes that require excessive computer time to provide quality 
results. 

In an attempt to bridge this gap between simplistic linear 
theory methods and complex Euler solvers, several 
methodologies such as the second-order potential analysis, 7 
hypersonic small disturbance theory, 3 and, more recently, 
nonconservative full potential methods 9 - 10 have been con- 
sidered by various investigators. The second-order theory, 7 in 
spite of the significant improvements reported, suffers from 
the lack of nonlinearity in resolving proper cross flow shocks 
arid sonic lines. Also, the singularities inherent in the for- 
mulation create difficulties in the numerical treatment of 
subsonic leading edges. The finite difference analysis of the 
hypersonic small disturbance theory 8 indicates that the 
solution procedure is as complex as that for the Euler 
equation and not particularly responsive to preliminary design 
level of effort. 

Recently, Grossman 9 and Grossman and Siclari 10 have 
computed supersonic fiowfields over conical and nonconicai 
cambered and twisted delta wings with remarkable success 
using the nonconservative full potential equation and a 
transonic relaxation method. However, their approach is 
made complicated by the use of global conformal mappings 
which apply only to certain classes of configurations. Also, 
the nonconservative form of the full potential equation is in 
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terms of second derivatives of the potential <t>, which, when a 
transformation is applied, generates a large number of first 
and second derivative transformation terms. 

The full potential method proposed in this paper is 
significantly different from that of Refs. 9 and 10. First of all, 
the method is based on the conservative form of the full 
potential equation, since for a shock capturing procedure to 
conserve mass across the shock wave," it is essential that the 
equation be cast in conservation form. 12 Second, the method 
can accommodate a numerical or analytical mapping 
procedure that is either orthogonal or nonorthogonal without 
complicating the form of the equation, in contrast to Refs. 9 
and 10. Third, the method is based on an approximate fac- 
torization implicit algorithm that can yield convergence much 
faster than the conventional successive line over-relaxation 
method. 10 Finally, the method is not an adaptation of a 
transonic code using type dependent operators, but a scheme 
specifically developed and tailored for supersonic marching 
problems using a density linearization concept and has no step 
size restrictions. 

To validate the present methodology, results are shown for 
a variety of conical and nonconicai geometries and are 
compared with Euler solutions and full potential results of 
Refs. 9 and 10. Results indicate that the method works just as 
fast and efficient for nonconicai flows as in the case of conical 
geometry treatment. Results also indicate that the method is 
very useful in computing very high-speed flows (A/„ -.2-6) for 
the moderate flow deflection angles (a-4-10 deg) where the 
neglect of entropy generation does not seriously distort the 
main features of the fiowfield. 

The present method can also handle more complicated 
geometries (realistic wing-body combinations) than the ones 
reported in the paper, but requires a suitable grid generation 
routine, especially near wing-body junction regions. In a 
subsequent paper, 13 results for nonconicai wing-body flows 
will be presented along with a formal method of charac- 
teristics treatment for cross flow signal propagation. 

II. Formulation 

The conservative form of the full potential equation in 
Cartesian coordinates x, y, z can be written as 

d(pu) d(pv) d(pw) 

dx dy dz ( 

where p is the density and u, v, w are the velocity components. 
They are calculated as the gradient of the potential <t>, 

u = <f> r ; v=4>/, w=<j>. (2) 
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IMPLICIT MARCHING SCHEME FOR SUPERSONIC FLOWS 


The density p is computed from the isentropic formula 

[ y — I ll/y-l 

1 -M 2 ,(u 2 + v 2 + w 2 -l)\ (3) 

If the density is normalized with respect to the freestream 
value, then the speed of sound a is given by 

a 2 =p T ' , /A/^ (4) 

where M m is the freestream Mach number. 

The objective of the paper is to solve for the scalar potential 
<t> from Eq. (1) subject to the surface tangency condition 
<t>„=0(n is normal to the body surface). Examining Eq. (1), it 
is very clear that <t> appears in a nonlinear form due to the 
presence of the density term inside the derivative. The ap- 
proach to be described here is a method that treats the density 
term in such a way that it produces the correct artificial 
viscosity needed for shock capturing and that enables one to 
solve for <t> with relative ease. 

In order to apply the surface tangency condition at the 
actual body location, a body-fitted coordinate transformation 
is essential. Introducing a body-fitted coordinate trans- 
formation, f=f(x, y.z). and $ = Z(x,y,z), Eq. 

(1) transforms to 

where U, V, and W are the contravariant velocity com- 
ponents. Introducing the following notation for convenience: 
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the contravariant velocities and density are given by 
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The Jacobian of the transformation J is represented by 


which indicates that for stability, the form of artificial 
viscosity be 


[' ■ - J 2 ) [AW V+m + I 

+ A| W\ u<t> n( + + W<t, m 1 ] (8) 


assuming that U, V, W are positive. What this implies is that if 
the flowfieid is hyperbolic (q>a), then solution can be ob- 
tained by marching along the hyperbolic flow direction s. 
Once the total velocity q becomes less than a, then marching 
along s is not possible. This is reflected in the fact that the 
effective artificial viscosity given by Eq. (8) is now negative. 

Now we will proceed with the numerical procedure for 
solving Eq. (5), and show the resemblance of the resulting 
artificial viscosity to that of Eq. (8). 


A. Treatment of in Eq. (5) 

Consider the direction f to be the marching direction. The 
condition to be satisfied for this to be true will become evident 
at the end of this analysis. Both the density p and the con- 
travariant velocity U are functions of the potential <t> and the 
transformation metrics, as represented in Eq. (6). In order to 
finite difference this f derivative quantity in terms of <t> only 
will require some linearization treatment of the density. This 
will be termed the “local density linearization” procedure. In 
the transonic formulation described by Holst, 13 the density is 
upwind biased and computed at the old level, while retaining 
central differencing for the ( U/J) term at the current level. 
Such an upwind density bias is shown to produce the right 
artificial viscosity in Ref. 14. Referring to Fig. 1, for a pure 
supersonic marching problem (say we want to march from the 
Ah plane to the /+lth plane), a transonic relaxation 
procedure 13 in the marching direction f is not appropriate 
because the solution <t> at the /+ 1th plane is not influenced by 
the /+2th plane. Hence, the following marching procedure is 
developed. 

Given the <t> information at all the previous planes i, /-l, 
i- 2 ,..., the problem is to compute <t> at the /+lth plane. 
Now, expand the unknown p=p(<t>) in terms of a neighboring 
known state denoted here by a subscript 0 (Ah plane in- 
formation would represent the neighboring known state for 
the/+ 1th plane). 


P=Po + 



A + 


(9) 


a<r.ii.t) 

d(x,y,z) 


'f, f, r ; 

Vr V, V; 

Jx €, . 


(7) 


Equation (5) is in terms of a general coordinate system 
(f.»7.S) and can accommodate any kind of mapping 
procedure, either analytical (conformal mapping) or 
numerical type. Any numerical marching procedure applied 
to Eq. (5) to simulate a supersonic flow should have a 
truncation error whose leading terms represent a correct 
artificial viscosity. This is essential to ensure marching 
numerical stability and to exclude the formation of expansion 
shocks which are unphysical and correspond to a decrease in 
entropy. The nature of the required artificial viscosity can be 
studied by an analysis 14 of the canonical form of Eq. (5), 
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Fig. I Implicit computational molecule. 
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where A<j> = (t>—<t> 0 and dp/d<t> can be shown to be a differential 
operator 16 



Substituting Eqs. (9) and (10) into the first term in Eq. (5), we 
get 


d(pUU) 

dS 




dS 




(11) 


Substituting for U in terms of <j> from Eq. 6 and rearranging 
Eq. (1 1) in terms of the potential difference A <t>, we get 
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where the speed of sound a 0 , the density p 0 , and the con- 
travariant velocities (J 0 , Vo. W 0 , represent information at the 
neighboring known plane. The f derivative term of Eq. (12) 
will now be one-sided differenced. Assuming U is positive. 
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An upwind differencing of the form Eq. (13) applied to Eq. 
(12) can be shown to produce a truncation term whose leading 
term is 



which will always represent a positive artificial viscosity as 
long as 

(J 2 

— >a 2 (15) 

a„ 

The preceding relation sets the condition for S to be the 
marching direction, for if U 2 /a n is less than the square of the 
local speed of sound, then the artificial viscosity becomes 
negative and a marching instability will occur. This also 
implies that the projection of the total velocity vector q in the 
direction normal to the const plane (t/,£ plane) is super- 
sonic. For example, in a spherical system, for the radial 
direction r to be the marching direction, the radial velocity q r 
must be supersonic. The similarity between the artificial 
viscosity term given by Eq. (14) and the first term appearing in 
Eq. (8) can be noted. When backward differenced, the terms 
in Eq. (12) will lead to a diagonally dominant tridiagonal set 
of equations for the unknown A <t> when coupled with the other 
two terms in Eq. (5). The mixed derivative terms like <t>-, and 
<5-e appearing in Eq. (12) will be upwind biased, depending on 
the sign of the coefficient multiplying them to preserve 
diagonal dominance and to provide the right artificial 
viscosity. 


d(pV/D 

B. Treatment of in Eq. (5) 

3i| 


This term will be written at the /+ 1th plane to make the 
resulting scheme fully implicit. 


The density term p in Eq. (16) cannot be represented at the 
/+ 1th plane since that would result in a very complicated 
nonlinear form for <t>. Hence, a density approximation is 
introduced by setting p = p* where p'=p 0 for conical flow 
treatment. In the case of nonconical flows, while advancing 
from /' to the < -H 1th plane, several iterations are performed 
within each cross flow plane (t),£) to refine the density p* to 
properly account for the axial geometry variation. This is 
done by initially setting p* to p 0 and then subsequently 
refining it by setting p* to the previous iterate value of p at the 
current plane. In many cases where the axial variation of the 
geometry is gradual (especially for smaller step size 
calculations) it was found that setting p* =p 0 even for non- 
conical flows produced very good results without having to 
refine the density subsequently. 

Writing Eq. (16) in terms of the potential difference A <t> 


3(pVU) 

. ( P*°z/ , 

fp‘a 

d-n 

"V J as)*' 

V J 



a 

— <t>. 

\ J 

a$ /, V / 

dv 


3tj 


a? 


A simple central differencing for the various terms in Eq. 
(17) will not be sufficient as that would not provide the 
desired artificial viscosity given by Eq. (8), required for shock 
capturing. To simulate an artificial viscosity of the form given 
by Eq. (8), the density will be upwind biased based on the 
previous work reported in Refs. 14-16. The density p* will be 
replaced by a modified density p* given by 


(p*);+«./t — U — v j+n.ic) (p*)y+'/,.t 

+ *4. 7 + w l ( / + «) (P*W* + (/ -«) (pV, + *u 108) 

where m = 0 when (V 0 ) j+ „, k >0 and m= + l when 
(Vo) /+•/,. k <0- When 6 is set to zero, first-order accurate 
density biasing is achieved while 6 = 2 gives second-order 
accuracy. The artificial viscosity coefficient , /l k is com- 
puted as follows: 

v i+'A.k = U~ (Po/^oJly+i.* (19) 

where s=0 for V J+l/i k >0 ands= 1 for Ky +//S it <0. 

Treatment of density as represented by Eqs. (18) and (19) 
would always produce a positive artificial viscosity as long as 
the local total velocity q 0 is supersonic. If that becomes 
subsonic, then the marching procedure would fail and the 
problem have to be treated as a transonic problem. 

The treatment of the (d/d^)[p\V/J\ term in Eq. (5) is very 
similar to the just described (d/d-q)[pY/J] term, except that 
the density biasing will be in the £ direction and will be based 
on the sign of fV. 

C. Implicit Factorization Algorithm 
Combining the various terms in Eqs. (12) and (17), and the 
terms arising from (d/d(){plV/J] will result in a fully implicit 
representation of Eq. (5) which cannot be solved without 
introducing an approximate factorization procedure. After 
some rearrangement of the terms, the factored implicit 
scheme becomes 
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This equation has the form 

L'L,(A<t>)=R (21) 

and it is implemented as follows: 

L s (A<t>)'=R I,(A4>) = (A<f>)‘ 4>=<f> 0 + A<i> (22) 

The various quantities appearing in Eq. (20) are given by 



and the right-hand side term R consists of various known 
quantities. 

The algorithm Eq. (22) requires only scalar tridiagonal 
inversions. Also, the scheme does not pose any restrictions on 
the direction of sweep that are present in the successive line 
over-relaxation method. 9 ' 10 


D. Freestream Truncation Errors 
To subtract out any numerical truncation error due to 
incomplete metric cancellation, 16 it is essential to add the 
terms (especially for a highly stretched nonorthogonal grid) 


a /p„c/„' 

\ + 3 

u 3 (^) 

3f \ J . 

' dr) \ / / 

> 3$ V J ) 


to the right-hand side of Eq. (20). 

E. Boundary Conditions 

In order to solve for A <t> from Eq. (20), boundary con- 
ditions will have to be prescribed at all four boundaries as 
shown in Fig. 2 at the current /+ 1th plane. While performing 
the L ( operator in Eq. (21), boundary conditions in terms of 
A<t>* will be required along the k= 2 and Ar=XMAX-l 
boundaries. For a pure angle-of-attack problem, k-2 and 
it = KM AX -1 can be considered as planes of symmetry 
across which all flow variables reflect. The quantity A<£’, even 
though it has no physical significance, can be safely set 


( A<t> ) T+ /j.XMAX ~(A4>)* + lj,KMAX-2 


U*)? +u/ = <A*)f +w (25) 

The L, operator would require boundary conditions along 
7 = 2 and j=JMAX in terms of Ad>. Since./' =2 is the body, the 
surface tangency condition 


V — a 2!<tl r ^ Q 22® t ^23^1 — 0 (26) 



I ■ 2 IS BODY 

I ■ 1. k » 1. k ■ KMAX ARE DUMMY ARRAYS 
I • JMAX IS FREE STREAM 

Fig. 2 Physical and computational plane. 


will be set at all points (t'+l,2,Ar). Along j=JMAX, 
freestream A <t> will be imposed. 

F. Grid System 

As shown in Fig. 2, the physical space ( x,y,z ) is trans- 
formed into a body-fitted (f,i 7 ,£) computational space. The 
transformation is performed numerically by using the elliptic 
grid generation techniques originally developed by Thompson 
et al., 17 and later modified by Steger and Sorenson, 18 and 
Middlecoff and Thomas. 19 The present full potential method 
does not require an orthogonal grid; however, the error in- 
troduced by the approximate factorization, Eq. (20), can be 
minimized if the grid is orthogonal in the cross flow plane 
(ij,|) . For conical flow calculations, the grid is generated only 
once, and, as the marching procedure continues, the grid is 
allowed to grow conically. For a general nonconical body, it 
would be necessary to construct the grid in every marching 
plane. 


III. Results 

Results are presented for both conical and nonconical 
supersonic flows. Comparisons are made with Euler 3,3 and 
full potential 9 - 10 results and experimental data. All the 
calculations were performed using a CDC 7600 machine. 

A. Conical Flows 

Besides validating the methodology, computation of 
conical flows is of interest for generating the initial data plane 
for nonconical calculations. For a conical geometry (radially 
invariant), the initial data plane with freestream conditions is 
chosen at some location f=f 0 (usually set at f=l). The 
solution is then marched along fusing Eq. (20) and boundary 
conditions. The conical flow calculation is assumed to have 
converged when the change in the root mean square density is 
less than 10 

Supersonic Leading-Edge Delta Wing 

Figure 3 shows the compression surface pressures for a 
supersonic leading-edge delta wing at M m = 6, angle of attack 
a = -8 deg, and leading-edge sweep A = 70 deg. The present 
full potential solution compares well with the Euler solution 19 
and experimental data. Also shown are the results from the 
first- and second-order linear theory. 7 Using a 30 x 40 grid in 
the (jj,£ ) plane, the present approach required 40 iterations to 
achieve convergence, and 12-15 s of computer time to produce 
the results shown in Fig. 3. 

It is interesting to note that in spite of the limitations of the 
full potential theory, even at a very high Mach number of 6, 
the comparison is in reasonable agreement with the Euler 



Fig. 3 Predicted compression surface pressure for a 70-deg sweep 
delta wing at = 6, a = - 8 deg. 


66 





V. SHANKAR 


AIAA JOURNAL 


t 

solution and is significantly better than the second-order 
theory. The discrepancy between the full potential and Euler 
results is mainly due to neglect of entropy generation in the 
present approach. 

Circular Cone and Ellipse 

Figure 4 shows the surface pressure distribution for a 
circular cone, half angle 7.5 deg, M m =3 at 15 deg angle of 
attack. At this angle of attack the cross flow Mach number 
becomes supersonic as the flow turns around the cone from 
the windward symmetry to the leeward symmetry. This cross 
flow supersonic region is terminated by the formation of an 
embedded shock on the cone surface. This is evident from the 
results of Fig. 4. Grid clustering near the cross flow shock was 
used both in the Euler calculation of Kutler, 3 and in the 
present method, to finely resolve the pressure jump. The 
present calculation required 25 s of computer time, using a 
30 x 60 grid in the (ij,{ ) plane. 

The liftoff of the vortical singularity on the leeward 
symmetry plane associated with the formation of the em- 
bedded shock is shown in Fig. 5. The location where the 
contravariant velocity V goes through zero on the leeward 
symmetry plane (1^=0) denotes the location of the vortical 
singularity. The behavior of the cross flow streamlines 
converging to the vortical singularity is also shown in Fig. 5. 

Figure 6 shows the full potential and Euler cross flow Mach 
number contours for the circular cone case. The presence of 
the embedded shock wave in both the results is very clear. The 
location of the vortical singularity liftoff is also shown in the 
figure. The Euler result is very oscillatory near the vortical 
singularity location while the present method predicts a 
smoother flowfield in the vicinity of the vortical singularity. 

The surface pressure distribution on an elliptic cone 
<><.= 18.39 deg, 5 C = 3. 17 deg at M m = 1 .97 and a= 10 deg is 
shown in Fig. 7. The results of the present study are compared 
with Euler calculations of Siclari, 5 full potential results of 
Grossman,’ and the linearized thin wing solution of Jones 
and Cohen. 20 The agreement between the various nonlinear 
methods is very good, including the position and strength of 
the embedded shock wave. 

Wing-Body Combination 

Figure 8 shows the numerically generated grid distribution 
in the cross flow plane (tj,£) of a conical wing-body com- 
bination. The design of this conically cambered delta wing to 
achieve shockless recompression is reported in Ref. 21. Figure 
9 shows the pressure distribution around this wing-body 
combination at M m = 2, and a = 7.81 deg. The leading-edge 
sweep A is 57 deg. The comparison of the results from the 
present method with the experimental data 21 is excellent. The 
calculation used a 15 x49 grid in the (ij,£) plane and required 
less than a minute of computer time. 




Fig. 5 Vortical singularity liftoff for a circular cone at =3, 
a= 15 deg, 9 C » 7.5 deg. 




Fig. 6 Comparison of cross 
flow Mach number contours 
for cone at angle of attack; 
Af„=3, ot= 15 deg, 9 C =1.S 
deg < VS = vortical singularity). 


PRESENT EULER 

FULL RESULT 

POTENTIAL (REF.3) 
RESULT 



Fig. 4 Surface-pressure distribution for cone at angle of attack; Fig. 7 Surface-pressure distribution on an elliptic cone; M m = 1.97, 
M m =■ 3, a ■ 15 deg, 9 C ■ 7.5 deg. 9 C = 18.39 deg, S c = 3. 17 deg, a = 10 deg. 
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Fig. 8 Computational grid around a wing-body combination. 



Fig. 9 Surface-pressure distribution for a conically cambered wing- 
body combination; Af, =■ 2, a «• 7.81 deg, A =■ 57 deg. 


B. Nonconical Flows 

Results are also presented for nonconical bodies of 
revolution and compared with experimental data. The initial 
data plane for the nonconical marching calculation is first 
obtained by performing a conical calculation over an assumed 
very small conical nose. This conical calculation usually takes 
20-30 iterations on a typical 30x30 grid in the (n,l) plane. 
The nonconical calculations did not exhibit any increase in 
computational time over the conical procedure. As mentioned 
earlier, for the applications considered here where the cross 
flow station does not vary substantially from the previous 
one, it was found that there was no need to iterate the solution 
at each cross flow plane ()?,£) to converge the density, and 
plottable accuracy was achieved by simply marching right 
along the body. However, if the body changes shape ap- 
preciably, the current implicit procedure might take 3-5 
iterations per cross flow plane to refine the solution. 

Reference 22 contains experimental data for several bodies 
of revolution at various Mach numbers and angles of attack. 
The shape chosen for comparison here is a circular arc- 
cylinder body. After the initial data plane was computed using 
a conical nose assumption, the current method typically used 
60 marching steps to reach the end of the body but the 
calculations are not subject to any step size restriction. A 
typical calculation required 40-45 s of computer time. Figure 
10 shows the circumferential surface pressure distribution at 
two different axial stations (x/l= 0.225 and 0.425) for 4 and 8 
deg angles of attack at M m = 2.3. The results from the present 
method are compared with experimental data, 22 showing very 
good agreement for the windward region with some 



Fig. 10 Circumferential pressure distribution for a circular-arc- 
cylinder body at M m » 2.3. 



Fig. 11 Surface-pressure distribution for a circular-arc-cyiinder body 
ala>8deg,M. >2J. 


discrepancy on the leeward side, possibly due to boundary- 
layer buildup. 

Figure 11 shows the surface pressure distribution in the 
axial direction along the windward and leeward plane of 
symmetry, at =2.3 and at = 8 deg. Again, results from the 
present method compare very well with the experimental 
data. 22 

IV. Conclusions 

An aerodynamic prediction technique based on the full 
potential equation in conservation form is developed for the 
treatment of supersonic flowfields. A local density 
linearization concept and a second-order accurate density 
biasing scheme are introduced in developing an implicit 
marching procedure. The method produces results that 
compare well with Euler solvers, and requires an order of 
magnitude less computer time and significantly less computer 
memory over existing Euler codes for the cases presented in 
the paper. In a subsequent paper, 15 results for more com- 
plicated nonconical wing-body flows are presented, along 
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with a formal theory for the characteristic signal propagation 
in the cross flow plane. 
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Abstract 

A nonlinear aerodynamic prediction technique 
based on the full potential equation in conservation 
form has been developed for the treatment of super- 
sonic flows. The method uses the theory of charac- 
teristic signal propagation to accurately simulate 
the flow structure, which includes shock waves and 
mixed elliptic-hyperbolic crossflow. An implicit 
approximate factorization scheme is employed to 
solve the finite-differenced equation. The neces- 
sary body-fitted grid system in every marching 
plane is generated numerically, using an elliptic 
grid solver. Results are shown for conical and 
nonconical wing-body combinations and compared with 
excerimental data and Euler calculations. The 
method demonstrates an enormous savings in execu- 
tion time and memory requirements over Euler 
methods. 

I. Introduction 

Nonlinear aerodynamic prediction techniques 
based on the Euler equation 1-3 and the full poten- 
tial equation 4-7 are steadily maturing into complex 
aerodynamic tools and becoming an attractive alter- 
nate approach to using the linearized panel 
methods. 3 Panel methods can handle very compli- 
cated geometries requiring minimal computer time 
to provide less accurate results, while the Euler 
solvers need expensive computer runs even for 
simple wing-body configurations. The full potential 
methods 5,7 are a substitute for the Euler methods 1 ’ 3 
to avoid the requirement of excessive computer time 
and memory allocation. While using a full potential 
method for supersonic flows, one should be aware of 
the isentropic limitations of the theory. As a 
general rule, the full potential theory is expected 
to perform well when the product of the Mach number 
and the characteristic flow deflection angle is less 
than one (M6 < 1) . 

The full potential method of Refs. 4-6 is based 
on the nonconservative form of the equation, while 
Ref. 7 and the present paper deal with the conser- 
vative form, to conserve mass across the shock. 4 ’ 1 " 
In order to properly treat the supersonic flow 
structure, which includes shock waves and mixed 
elliptic-hyperbolic crossflow, the present method 
uses the theory of characteristic signal propaga- 
tion based on the eigenvalue system of the full 
potential equation. An approximate factorization 
implicit scheme that includes a density biasing 
procedure in the crossflow plane activated by the 
eigenvalue system is incorporated to accelerate the 


♦Member Technical Staff, Associate Fellow AIAA 
♦•Professor, Department of Mathematics 


computational efficiency. The implicit scheme does 
not pose any restrictions on the direction of 
sweep that are present in the successive line 
overrelaxation method (SLOR). 4- '' 

The full potential as well as Euler methods 
require the application of boundary conditions at 
the actual body surface location. This, in general, 
necessitates the use of a body-fitted coordinate 
system. In the present method, the equation is cast 
in a most general arbitrary coordinate system and 
the appropriate body-fitted grid is generated 
numerically, employing an elliptic grid solver. 11 

The paper presents various results for conical 
and nonconical wing-body configurations and com- 
parison is made with experimental data and Euler 
solution. The effect of the density biasing based 
on the characteristic signal propagation is demon- 
strated in terms of a sharper pressure profile 
across the shock wave. References 5 and 6 present 
excellent results at low supersonic Mach numbers, 
while Reference 7 and the present paper demonstrate 
the capability of the conservative full potential 
approach in handling even very high Mach number 
flows (Moo - 4-6, a - 0-8°). All the calculations 
reported in this paper were performed using the 
CDC 7600 machine and clearly demonstrated an order 
of magnitude or more reduction in computer time 
over Euler methods. A typical nonconical wing-body 
calculation takes less than 2 minutes of execution 
time to produce results comparable with experimental 
data. 

II. Formulation 

The conservative full potential equation cast 
in an arbitrary coordinate system defined by, • 

C = ?(x,y,z), n = n(x,y,z), and Z = ?(x,y,z), takes 
the form 


B) e * 

H)n 

* (»S) f ■ 0 • (“ 

where U, V, and W are the 
components. Introducing 
for convenience 

contravariant velocity 
the following notation 

U - Uj , 

V = u 2 

, W = U 3 

rH 

X 

II 

X 

y * * 2 

, Z = Xg 

5 ■ Xj , 

n * X 2 

, z = *3 


the contravariant velocities and density are 
given by 
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■( 2 ) 


P - [l - + v* n + w$ 5 - 1^1/(Y-1) 


The Jacobian of the transformation J is represented 
by 


3(c,n,S) 

J = = 

3 (x,y,z) 



' (3) 


Equation ( 1 ) is in terms of a general coordinate 
system (q,n,q) and can accommodate any kind of map- 
ping procedure, either analytical (conformal map- 
ping) or numerical type. Use of Eq. (1) to simu- 
late the supersonic flow by marching in the q 
direction, first requires the establishment that 
the equation is indeed hyperbolic with respect to 
the marching direction. The nature of Eq. ( 1 ) can 
be analyzed by studying the eigenvalue system of 
Eq. (1). Combining the irrotationality condition 
in the (q,n) and (q,q) plane and Eq. (1), one can 
write the following matrix equation. 


Aq c + Bq n + Cq ? = 0 


(4) 


where 


A = 


j <»»>, i (»“), j (»»), 


B 





-10 0 


0 0 



The subscripts in Eq. (4) denote differentiation 
with respect to that variable. 

In order for Eq. (4) or Eq. (1) to be hyper- 
bolic in q direction the following two conditions 
must be satisfied. 

( 1 ) A " 1 must exist 

(2) A^l real linear combinations of A _1 B and 
AC must have real eigenvalues (character- 
istics). This implies A -1 (aB + SC) must 
have real eigenvalues for all combinations 
of a and B satisfying a 2 + B 2 = 1 . 

When the above two conditions are applied to Eq. (4), 
the following criteria is obtained for q to be the 
marching direction. 



Where the transformation metric an is defined in 
Eq. (2) and a is the local speed of sound. Equation 
(5) is the most general form. For example, in a 
spherical system (r,e,$), for the radial direction 
r to be the marching direction, according to Eq. 

(5), the radial velocity q f must be supersonic. In 
a Cartesian system (x,y,z), for x to be the march- 
ing direction, the velocity u has to be supersonic. 
For convenience, the derivation of Eq. (5) for a 
Cartesian system is described in Appendix A, and 
the derivation for an arbitrary coordinate system 
(q.n.S) has been derived in a similar manner. 

Thus far, the condition for q to be a marching 
direction has been identified from the character- 
istic theory. Thes means the (n,C) plane will be 
treated as a marching plane, which will be defined 
from now on in this paper as the crossflow plane 
(the real crossflow is the projection of the veloc- 
ity vector on a unit sphere with center at the 
origin). Even though, the flow is supersonic in 
the marching direction (i.e. , hyperbolic type), 
the behavior of the flow structure in the cross- 
flow plane (n.S)'can be a mixed elliptic-hyper- 
bolic type. Depending on the nature of the flow at 
a crossflow plane grid point (whether elliptic, 
parabolic or hyperbolic), the n and £ derivative 
terms in Eq. (1) will be appropriately modeled. 
Again, the theory of characteristics will dictate 
how the signals are propagated in the crossflow 
plane. 


C 


1 (pw)* j (pW) 3 (pW)„ 

C n 5 

0 0 0 
-10 0 


A. Crossflow Signal Propagation 

The nature of the flow in the (n,C) plane can 
be analyzed by separately studying the eigenvalues 
of A _1 B and A _ 1 C. The eigenvalue character o'f A _1 B 
will determine the n- derivative treatment and 
similarly A _1 C for the I derivative. For illustra- 
tion, only the study of A _1 B is shown here, and 
A _1 C follows the same procedure. 
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The eigenvalues X of A _1 B are obtained by solv- 
ing (A _1 B - XI) = 0. Since A -1 is assumed to exist 
(condition (1) before Eq. (5)), the following is 
true. 




-1 

0 


0 


° ( 6 ) 


Solving for X, 


1,2 


(oU) ♦ (pV) 
n 




(PU) Mpv) 


- <(PU)* <PV>* 


2 { pU ) . 


(7) 


hyperbolic, appropriate finite difference models for 
the n derivative term in Eq. (1) are chosen. This 
will be described later in this paper. 



CASH 1. ELLIPTIC 

CROSSFLOW 



CASE 2. HYPERBOLIC 
CROSSFLOW 
POSITIVE V 


Now, analyzing Xj and X2 the following combinations 
are possible. 

(1) X^ is positive Xj^ is negative 

X 2 is negative X 2 is positive 

(2) Xj and x 2 both positive 

(3) Xj and X2 both negative 

(4) Xj or X 2 Is zero 

These possible combinations are schematically shown 
in Fig. (1). Each one of these combinations describe 
a different feature of the flow in the crossflow 
direction n. Referring to Fig. (1) diagram the 
following descriptions are made. 



CASE 3. HYPERBOLIC 
CROSSFLOW 
NEGATIVE V 



CASE 4. PARABOLIC 
CROSSFLOW 


Fig. 1. Eigenvalue structure in (c.n) plane. 


Case 1 

Here, one eigenvalue is positive and one negative 
This implies an elliptic type crossflow because the 
characteristic signals are brought into the point P 
from both the positive and negative direction of n. 

Case 2 

Here, both the characteristics are positive, 
which means the characteristic signals propagate 
into the point P only from below and anything hap- 
pening above the point P doesn’t influence that 
point. This describes a hyperbolic type crossflow 
point with a positive contravariant velocity V. 

Case 3 


One can readily see from Eq. (7), that one of 
the eigenvalues go to zero when 

(pV) =0 (8) 

% 

From the definition of p and V from Eq. (2), one 
can write 

= p | a 22 "a 2 ! {9) 

v 2 

Thus, when a 2 2 = occurs, the method will antici- 
pate a switch in the character of the crossflow, 
and realize the onset of the formation of a super- 
critical crossflow. 


Here, the characteristic signals propagate 
from above into the point P and similar to case 2, 
this describes a hyperbolic type crossflow with a 
negative contravariant velocity V. 

Case 4 


Besides providing a valuable information regard- 
ing the type of crossflow, the eigenvalues \, and 
X 2 of A _1 B and similarly \. and X. of A" l C can also 
be used to determine the Courant number in n and 
direction for a given step size Ac. 


Here, one of the eigenvalues is zero and 
describes a parabolic type crossflow. This will 
represent the crossflow sonic line. 

The transition from an elliptic to a hyperbolic 
crossflow type takes place through a parabolic 
point, which is indicated by one of the eigenvalues 
going to zero. Thus, by monitoring the eigenvalues 
Xi and X2 one can precisely model the crossflow 
plane terms. Depending on whether it is elliptic or 


(xmax) 7^- = Courant Number in n direction. 

' 1 n An 

(xmax) ^ ||- = Courant Number in I direction. 

Where (Xmax)n and (\max)^ define the maximum of 
(Xj ,X 2 ) and (X 3 ,X 4 ) respectively. 
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B. Treatment of (pj);. In Eg- (1) 

The direction c has been identified to be the 
hyperbolic marching direction satisfying the con- 
dition given by Eq. (5). Referring to Fig. (2), 
this derivative term will be backward differenced 
as, 


/11X ( a i * ■ h i)|W) <+1 * W), }" "’lit®, * 

‘ a 1 « 1 - Bbjjacj + ac 2 ) 


1-1 


( 10 ) 


where 


a l = ( A? l + A? 2 ) J 
b l " ( Ac l) 2 


9 = 0 first order accurate 


KNOWN DATA 

PLANE ^CURRENT 

OATA PLANE 
WHERE 0 IS UPDATED 



i+U.k + l 


i + l.j-l.k 


Fig. 2. Implicit computational molecule 


= 1 second order accurate 

Given the velocity potential <j> Information at 
all previous planes 1, 1-1, 1-2, .... the problem 
is to compute $ at the current plane 1+1. Equation 
(10) Involves both the density and contravariant 
velocity at the (i+1) plane, and both are functions 
of <£ (Eq. 1). In order to write Eq. (10) In terms 
of <p only will require a local linearization proce- 
dure. This Is done as follows: 


C. Treatment of the Crossflow Term /pjj in Eq. (1) 


Similar to the treatment of the (pU) term, the 
(pV).j + j will also be linearized as, 

(pV) 1+1 - (pV) 1 + [(PV)J. AO + ... (13) 


(pU) i+ i 4 (pu) 1 + JtpU^AO + ... (11) 

where 

(pu>* - + 



AO a ^1+1 " 

Substituting for p, and II, into Eq. (11), and group- 
ing various terms, ? 



The above locally linearized equation involves only 
A$ as the unknown to be solved for. To maintain the 
conservative differencing, both (pU)-)+i and (pU)-j 
appearing in the first term of Eq. (10) will be 
linearized. That is, (pU).,- will be linearized 
about (i-1) plane values. The upwind differencing 
of the q derivative term as shown in Eq. (10) will 
produce a truncation term whose leading term is 

)l — ( U 2 $__, At'. This will always represent 

( U 2 ) 

a positive artificial viscosity as long as the 
marching condition dictated by the characteristic 
theory, Eq. (5), is satisfied. 


The above linearized expression for (pV)i + i will be 
plugged inside the n derivative term of Eq. (1). 

It involves only AO as the unknown variable. The 


finite difference model for the 


(4 


term will be 


dictated by the theory of characteristic signal 
propagation as described in Section A of this paper. 
When the eigenvalues of A' 2 B represent case 1 in 
Fig. 1 (one positive and one negative eigenvalue 
representing an elliptic type), then all the terms 

in /pr\ will be central differenced. For this 


“■ ( a « ' s) 


is positive, and central differ- 


encing of the term along with the backward 

differencing of the |pjj term as in Eq. (10) will 


preserve the diagonal dominance. For cases 2 and 
3 of Fig. 1, the crossflow behaves like hyperbolic 


type, and ^a 22 - jzj is negative. Then, central 

differencing of the terms in Eq. (13) is inappro- 
priate, as it will destroy the diagonal dominance, 
and also will not provide the necessary artificial 
viscosity to avoid the formation of expansion 


shocks. Thus, when Xj and 


X 2 are both positive or 
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both negative (hyperbolic type), the terms in 
should be upwind differenced depending on 

the direction of V. However, such an upwind dif- 
ferencing in the n direction will not give rise to 
a tridiagonal system, and in general, the overall 
system will be pentadiagonal in nature. In order 
to preserve the tridiagonal nature of the implicit 
scheme, rather than upwind differencing the $ 
derivatives, the density biasing concept 7 ’ 12 is 
implemented when the crossflow is hyperbolic. 

The procedure is as follows: 


be the value at the i-th plane and then subse- 
quently iterated to convergence by setting p* to 
the previous iterated value of p at the current 
i+1 plane. 

A similar density biasing procedure is imple- 
mented for the term in Eq. (1). 

Activating the density biasing based on the eigen- 
value structure of A _1 B and A^C has proved to be 
very efficient in predicting sharp shock profiles 
The same concept can also be employed for tran- 
sonic applications. 


(pj) n * h t j ( a 21*e + a 22*n + a 23*s ) \ ( 14 } 

Here, the density p has been replaced by p defined 
to be (referring to Fig. 2) 

p 1+l,j+i,k (* ' ‘Vi.k) p *j+i,k + I'j+i.k 

jU+e) P*j +2m ,k +(1 ' e) p *j-l+2m,k} (15) 


where m=0 when 


(v n \ >0, = +1 when 
V °/j+i,k 

When 9 is set to zero, first 


C°W ‘ 

order accurate density biasing is achieved while 
9=2 gives second order accuracy. The artificial 
viscosity coefficient v... . is computed as 
follows: J *’ 


v 


j+i,k 



,j+s,k 


(16) 


where, q is the local total velocity, a is the 
speed of sound and 

s=0 for V j+Jjk > 0 

=1 for V j+i,k ' 0 

u=0 for - T-) > 0 (elliptic 

’ cc a /i,j+i,k crossflow) 

=1 for (a„ - ^-) <0 (hyperbolic 

\ a 'i,j+i,k crossflow) 


Combining the various terms of Eq. (1) as 
represented by Eqs. (10), (14) and (15) together 

with the terms arising from will result in 

a fully implicit model. This is solved using an 
approximate factorization implicit scheme. The 
details are given in Ref. 7. The boundary condi- 
tion at the body surface is the flow tangency 
condition, which is easily applied by setting the 
contravariant velocity V to zero at all body 
points. 


III. Grid System 

The transformation of the physical space 
(x,y,z) to a body-fitted computational space 
(c»n,S) is performed numerically by using the 
elliptic grid generation technique of Ref. 11. 

The body geometry at every marching plane is 
prescribed along with a suitable outer boundary 
where free stream conditions are imposed. Since 
the equation is cast in a general coordinate 
system, the marching plane (constant c) can either 
be a constant x plane or a spherical (constant r) 
plane as long as the marching criteria (Eq. 5) is 
satisfied. Given the geometry shape and the pre- 
scribed outer boundary, the following set of 
elliptic equations are solved to generate the 
interior grid. 

5 yy + 5 zz = P(C ’ n) 

(17) 

n yy + n zz = Q(C ’ n) 

The forcing terms P and Q are properly chosen 
to achieve two main desirable features: (1) to 

cluster grid points to a boundary, and (2) to 
force grid lines to intersect the boundary in a 
nearly orthogonal fashion. 


Thus, the density biasing is switched off when the 
eigenvalues \\ and A? exhibit an elliptic cross- 
flow. All the if> derivative terms are central 
differenced in Eq. (14). Treatment of the density 
as represented by Eqs. (15) and (16) would always 
produce a positive artifical viscosity when the 
crossflow is hyperbolic. The local total velocity 
is always assumed to be greater than the speed of 
sound, otherwise the marching procedure would fail. 

In Eq. (15), the evaluation of p* depends on 
whether the flow is conical or nonconical. For 
conical flows all p* quantities are evaluated at 
the i-th plane. For nonconical flows, at each 
nonconical marching plane, initially p* is set to 


Once the grid is generated, all the metric 
terms a^- in Eq. (2) and the Jacobian J in Eq. 

(3) are computed by numerical differentiation. To 
subtract out any numerical truncation error in the 
free stream due to incomplete metric cancellation, 13 
it is essential to add the terms (especially for a 
highly stretched nonorthogonal grid). 



to the right hand side of the finite differenced 
model of Eq. (1). 
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IV. Results 


A series of calculations were performed for 
conical and. nonconical geometries at various Mach 
numbers (M ^2 to 6) and angles of attack 
(a 'v 0 to 10°), to validate the full potential 
characteristic switch methodology, and assess the 
feasibility of using numerical grid solvers for 
complex configurations. The results from this 
study are compared with experimental data and 
Euler simulation. 

The generality of the formulation allows one 
to choose any c as the marching direction, provided 
the condition given by Eq. (5) is satisfied. 

Thus, depending on the geometry definition and 
the flowfield character one could choose either a 
constant x-plane marching or constant r-plane 
spherical marching. 

The effect of on, off density biasing based on 
characteristic signals in the crossflow plane 
(n,5), described by Eq. (16), is demonstrated in 
terms of the crossflow Mach number (M(;p) distribu- 
tion in the shock region in Fig. 3. When density 
biasing is applied everywhere 7 , including at 
elliptic crossflow points, it introduces unnecessary 
artificial viscosity and tends to smear the dis- 
continuities like shocks in the flow field. This 
is seen by the dash line crossflow Mach number 
distribution across the bow shock and across the 
embedded shock on a cone surface in Fig. 3. When 
the density biasing Is switched off at crossflow 
.elliptic points, the shocks appear as a sharper 
discontinuity (usually within two mesh intervals) 
as shown by the solid line distribution in Fig. 3. 
All the calculations to be presented here, were 
achieved using the second order accurate implicit 
scheme (8=1 in Eq. 10), with on,off density biasing 
activator u in Eq. 16. 

Figure 4 shows the grid arrangement in the 
marching plane for a conically cambered wing-body 
combination. The elliptic grid solver with orthog- 
onality constraints near the surface, required 40 
to 60 iterations to converge to within 10* 8 error in 
the residual. Figure 5 shows the pressure distribu- 
tion at M I = 2 and angles of attack of 7.81° and 
10.82°. The leading edge sweep is moderate (57°), 
and spherical plane marching is implemented (instead 
of x-plane marching) to avoid low supersonic Mach 
number components along x-direction near the leading 
edge. The results are compared with experimental 
data given in Ref. 14. The comparison is excellent. 
The marching step size A? is chosen by monitoring 
the eigenvalues and setting the Courant number to 
about 20. The numerical formulation, being a con- 
servative form, predicts a stronger. crossflow recom- 
pression on the leeward side than the ones seen in 
experiments. On a 20 x 49 (n,£) grid, the method 
requires about a minute of CDC 7600 time. The con- 
ical flow field is assumed to have converged when 
the change in root mean square density in the march- 
ing plane is reduced to less than 10 -5 . 

Figure 6 shows the surface pressure distribu- 
tion on a flat conical wing-body (that is not 
designed to weaken the crossflow shock formation) 
at two different angles of attack (1.72° and 5.71°) 
and Mach number of 2. The experimental data and 
the numerical prediction are in excellent agree- 
ment and clearly indicate the presence of an 
embedded crossflow shock. 


Even though the full potential theory is 
restricted by the isentropic assumption, one will 
be surprised to find out that the theory can be 
effectively utilized to predict even very high Mach 
number flows as long as M6 is less than or of the 
order of one (M6 £ 1). This is demonstrated in 
Fig. 7, which .shows the results for a Sears-Haack 
body at M„ = 6 and different angles of attach 
(0°, 4°, 8°). The numerical prediction is com- 
pared with an unpublished NASA Langley data, and 
the agreement is excellent. Constant x-plane 
marching is implemented for this configuration. 



ON.OFF DENSITY 

BIASING ACTIVATOR 
BASED ON 
CHARACTERISTICS 


DENSITY BIASING 




Fig. 3. Effect of density biasing activator on 
the crossflow Mach number distribution 
in the shock region 
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Fig. 4. Grid arrangement in the marching plane for 
a conically cambered wing-body combination 
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Fig. 5. Surface Dressure distribution on a 

conically cambered wing-body combination 

Moo = 2 


marching plane. The grid at each marching plane is 
generated using the elliptic grid solver. Figures 
9a to 9d show series of results at x/J. of 0.3, 

0.5, 0.65 and 0.8 respectively. The full potential 
results are compared with the experimental data and 
Euler simulation in Ref. 15. Figure 9a, which 



X 


Fig. 6. Surface pressure distribution on a flat 
conical wing-body combination Mo, = 2 




Figure 8 shows a schematic of a symmetric arrow 
wing-body configuration. The actual geometry shape 
is prescribed analytically as detailed in Ref. 15. 
Series of computer runs were made for this configu- 
ration at different Mach numbers and angles of 
attach and some results are presented here. First, 
an initial data plane near the nose region of the 
configuration is established by assuming a conical 
nose shape. The nonconical marching is then 
initiated. At each nonconical marching plane, the 
density is iterated to convergence (p* in Eq. 15 
usually takes 2 to 3 cycles, to converge to 10 " 5 
error tolerance) before proceeding to the next 


C P 





o 


Fig. 7. Pressure distribution on Sears-Haack body 
at Mo, = 6. Windward plane of synmetry. 
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shows results for x/l of 0.3, clearly demonstrates 
the accuracy of the full potential simulation. It 
is surprising to see that the present full potential 
method compares with the experimental data even 
better than the Euler calculation even at a high 
Mach number of 4.63. Similar excellent full poten- 
tial results are shown in Fig. 9b for x/l of 0.5 
and compared with data from Ref. 15. The striking 
full potential results are shown in Fig. 9c, where 
the unphysical oscillations experienced by the 
Euler simulation at H* = 2.36 near the wing-body 
junction area are not seen in the present method 
and comparison with experimental data is more 
dramatic. Figure 9d shows the pressure distribu- 
tion at x/l of 0.8, where the wing is separated 
from the body. The wake is simulated by assuming a 
planar shape, and imposing pressure equality (in the 
present method, it will be density equality due to 
full potential formulation) across the cut. Again, 
the full potential results are in good agreement 
with the Euler solution and experimental data. 

Figure 10 shows an angle of attack case, 

M = 4.63, a ■ 3° for the same symmetric arrow 
wing-body configuration. The results are compared 
with the data of Ref. 15 at x/l of 0.65, and the 
agreement is good even near the wing body function 
region. x 

A typical arrow wing-body calculation using a 
20 x 49 grid in the (n,E) plane and choosing a 
marching, step size Courant number of 3, required 2 
minutes of CPU time for the entire calculation, 
which includes the numerical grid generation at each 
plane and the conical initial data plane. This rep- 
resents an enormous savings in computer execution 
cost over other nonlinear methods, especially Euler 
solvers. 



X 

Fig. 8. Top and side views of a typical arrow 
wing-body configuration 





Fig. 9a. Grid arrangement and surface pressure 
distribution for a symmetric arrow 
wing at x/l = 0.3 
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Fig. 9b. Arrow wing pressure distribution at 
x/i = 0.5 


Fig. 9c. Arrow wing pressure distribution at 
x/i = 0.65 
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Fig. 9d. Arrow wing pressure distribution at 
x/i. = 0.8 
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Fig. 10. Angle of attack solution for the 
arrow wing at x/t = 0.65 
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V. Conclusions 


1) A" 1 exist 


A nonlinear full potential aerodynamic pre- 
diction capability based on a sound mathematical 
theory of characteristic signal propagation has 
been developed. The method uses a general body- 
fitted coordinated system and numerical mapping 
techniques. The on, off density biasing activator 
in the crossflow plane has proved to be very effec- 
tive in capturing sharp shock profiles. Results for 
conical and nonconical flows at various Mach numbers 
and angles of attach are shown to be in excellent 
agreement with experimental data and Euler results. 
The enormous savings in computational cost exhibited 
by the present approach, makes it a very promising 
substitute for the less accurate linearized panel 
methods, and expensive Euler solvers for use as a 
preliminary design tool. The future work will 
involve automatic grid generation for wing-body- 
nacelle-canard configurations and better wake 
treatment. 


2) A -1 (aB + 8C) must have real eigenvalues 
for all a and 8 satisfying a 1 + 8 2 = 1 

Since A' 1 is assumed to exist, the eigenvalues 
of A -1 (aB + 6C) can be obtained by solving 


(aB + 6C - XA) = 0 (A2) 

Substituting for A, B and C from Eq. (Al) into Eq. 
(A2), the roots of the equation are obtained. 



- ^va+wBl^fvaW-)!-^) |l- 



(va+w8) 2 

a 2 


(A3) 


Appendix A 

Derivation of the marching condition (Eg. 5) for a 
Cartesian system 


Eq. (A3) will have real values as long as the square 
root term is real. This implies, the quantity 
inside the square root must be positive. Simplify- 
ing the quantity inside the square root, the condi- 
tion becomes 


The Cartesian system analog of Eq. (4) is given 
by u 2 + (va + w8) 2 Q 

Af x + Bf + Cf 2 = 0 (Al) a " 


(A4) 


where, 

[7pu) u (pu) v (pu) w 


Let 


a = cos 6 
8 = sin 9 


=j» a 2 + 6 2 = 1 


0 0 1 
(pv) u (pv) v (pv) w 


B = 


-1 


0 0 


0 0 0 


v = q cos 9 
w = cf sin "9 

where , 



tan 9 = w/v. 

Substituting these into Eq. (A4) and simplifying 
results in 


C = 


(pw) u (pw) v (pw) w 


-1 


u 2 + q 2 cos 2 (6 - 9) 


- 1 >0 


(A5) 


Since this condition must hold for all combinations 
of 9 and ¥, Eq. (A5) implies (for ¥ - 9 = k/Z) 


f = 




for x to be the 
’ marching direction 


(A6) 


Equation (Al) is hyperbolic with respect to Eq. (A6) is a special case of Eq. (5) in Section II. 

x-direction if 
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geometries to evaluate the capability of the approximate equations of motion 
considered. Results from the computations indicate good agreement with higher order 
solutions and experimental results for a variety of wing, body, and wing-body 
shapes for ■values of the hypersonic similarity parameter MS approaching one. Case 
computational times of a minute were achieved for practical aircraft arrangements. 
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